% This code computes the minus of sum of log likelihood
function likeli = rstability_likelihood(rs_para, s2_mkt_mm, delta_p_mm, ...
    in_regression_input, reg_after, fc_value, N_st, mkt_input)

% parameters being estimated
if fc_value==1
    cost2_0     = rs_para(1); % scalar
    mu_cost2    = rs_para(2);
    sigma_cost2 = rs_para(3);
elseif fc_value==2
    cost2_0_pre   = rs_para(1); % scalar
    cost2_0_post  = rs_para(2); % scalar
    mu_cost2    = rs_para(3);
    sigma_cost2 = rs_para(4);
elseif fc_value==N_st
    cost2_0_mkt = rs_para(1:end-2);
    mu_cost2    = rs_para(end-1);
    sigma_cost2 = rs_para(end);
end

if fc_value==1
    cost2_0_mat = cost2_0*ones(size(s2_mkt_mm)); % NNx KM
elseif fc_value==2
    cost2_0_vec = (reg_after==0)*cost2_0_pre + (reg_after==1)*cost2_0_post;
    cost2_0_mat = repmat(cost2_0_vec, 1, size(s2_mkt_mm,2)); % NNx KM

elseif fc_value==N_st
    cost2_0_vec = zeros(size(s2_mkt_mm,1), 1); % NN x 1
    for ii=1:size(s2_mkt_mm,1)
        mm=mkt_input(ii);
        cost2_0_vec(ii)=cost2_0_mkt(mm);
    end
    cost2_0_mat = repmat(cost2_0_vec, 1, size(s2_mkt_mm,2)); % NNx KM
end

% likelihood when delta_p=0, pr0: N.ist x K*M
pr0 = 1-logncdf((s2_mkt_mm.^2)./(2*cost2_0_mat), mu_cost2, sigma_cost2);

% likelihood when delta_p>0, pr1: N.ist x K*M
pr1 = logncdf((s2_mkt_mm.^2)./(2*cost2_0_mat), mu_cost2, sigma_cost2)... % this is (1-pr0)
    .*lognpdf(s2_mkt_mm./delta_p_mm, mu_cost2, sigma_cost2);

% sum of log likelihood
likeli = sum(log(pr0(delta_p_mm==0))) + sum(log(pr1(delta_p_mm>0)));

% take minus
likeli =  - likeli/size(s2_mkt_mm,1);
